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Abstract. - We study non-equilibrium velocity fluctuations in a model for the sedimentation 
of non-Brownian particles experiencing long-range hydrodynamic interactions. The complex 
behavior of these fluctuations, the outcome of the collective dynamics of the particles, exhibits 
many of the features observed in sedimentation experiments. In addition, our model predicts 
a final relaxation to an anisotropic (hydrodynamic) diffusive state that could be observed in 
experiments performed over longer time ranges. 



Despite the study of sedimenting suspensions has a long and well-deserved history for their 
ubiquitous nature and applications [[l] , attention to the non-equilibrium density and velocity 
fluctuations in these systems has only been paid lately. In particular, the nature of non- 
equilibrium fluctuations in the sedimentation process has been a subject of a long controversy. 
While theoretical arguments ^ and extensive computer simulations m suggested that velocity 
fluctuations should diverge with the system size, the available experimental results ||, || , and 
the theoretical analysis in Ref. [0], found no evidence for such divergences. These apparently 
contradictory observations may have found a reasonable interpretation after the experimental 
evidence in Ref. (7) , and the theoretical study by Levine et al. [|| . 

Another striking piece in the puzzle of sedimentation was recently added by the experimen- 
tal work of Rouyer et al. ]£[. In their experiment, the authors analyzed the trajectories and 
velocities of non-Brownian [ H| particles sedimenting in a quasi-two dimensional (2d) fluidized 
bed, and showed the intrinsic non-Gaussian nature of velocity fluctuations. The main conclu- 
sions of this work are the non-Gaussian form of the probability density functions (PDF's) of 
the velocity fluctuations; the anisotropic character of the particle trajectories (diffusive along 
the horizontal direction and superdiffusive along the vertical one); and the presence of very 
long-range correlations in the velocity fluctuations along the gravity direction. New evidence 



along some of these lines is also provided in a recent paper by Cowan et al. 1 1 1 



The results of Refs. ||7|, |9|, |llj] pose new questions regarding the process of sedimentation, 
which have not been addressed by previous theoretical approaches. Our purpose in this Letter 
is to tackle these questions from the point of view of the particle's dynamics to ascertain 
the chief physical mechanisms underlying such fluctuation phenomena. In order to do so, 
we propose a model of sedimentation in which particles experience long-range hydrodynamic 
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interactions. We start from the solution of the linear Navier-Stokes equation for the suspension 
in an unbounded incompressible fluid . We consider a system of N particles obeying a system 
of coupled differential equations which we solve numerically. In the solution of the equations, 
we keep track of both positions and velocities of the particles, and compute several relevant 
statistical properties. In our model, we observe most of the experimental features reported in 
Refs. J?], |j| Q, namely, slow and fast particles, and swirls and channels in the velocity field, 
which, in overall, yield non-Gaussian velocity distributions and a slow time relaxation of the 
velocity autocorrelations. In addition, our model predicts a final relaxation to an anisotropic 
(hydrodynamic) diffusive state, not observed in Ref. M. 

The velocity of a particle n in a dilute suspension is given by the expression U„ = 
^2— H„ m • F m , where the sum is carried out over all the particles m in the suspension, H„ m is 
the mobility tensor, and F m , the external force acting on each particle, is gravity g, oriented 
along the positive z-axis jlj . The simplest form of the tensor H corresponds to dilute suspen- 
sions of point-like particles. In this case, the solution of the stationary Navier-Stokes equation 
in an unbounded medium yields the so-called Oseen tensor Ho(r) = (I + r ® r/r 2 )/8irr}r 
where I is the identity matrix, the operator <S> stands for the tensorial product, and r\ is the 
fluid viscosity. 

We study a suspension of monodisperse non-Brownian particles (for which incrtial effects 
are irrelevant in a viscous fluid) at very low concentrations, where the point-particle assump- 
tion is indeed a good approximation. Initially, particles are placed on the same vertical xz 
plane. The form of the mobility matrix in the Oseen approximation ensures that particles 
in such a configuration will never leave that plane. Simulations are performed on a system 
of N particles in a square cell of size L (corresponding to a concentration c = N/L 2 ). Pe- 
riodic boundary conditions (PBC's) are imposed in all directions (including the y direction 
perpendicular to the initial plane) in order to guarantee the uniformity of the suspension 
To avoid the discontinuities arising from truncating long-range hydrodynamic inter- 
actions, imposing PBC's amounts to considering Oseen interactions with an infinite set of 
images of the original system |i~4]]. In this way, the velocity of each particle is written as 
U ra = ^2 m Y2d Ho( r 7im + d) • g, where the index m runs through all the particles inside a cell 
of volume V, r nm indicates the relative position of a pair of particles within that cell, and d 
runs through the positions of the images of m in an infinite number of cell replicas along the 
x, y, and z axes. 

Imposing PBC's along the y axis is mathematically equivalent to imposing slip boundary 
conditions to the fluid velocity field on effective walls parallel to the sedimentation plane, and 
located at distances ±L y /2, where L y < L. Sedimentation experiments are usually carried 
out within a thin fluid slab confined by parallel glass plates. A realistic modelization of this 
process should thus include wall effects by imposing no-slip boundary conditions on the walls. 
By doing so, hydrodynamic interactions become exponentially screened for length scales larger 
than the slab thickness L y . One then expects that exponentially damped interactions intro- 
duce an external characteristic length into the problem, L y , which will govern the dynamics of 
the system, a fact that has not been pointed out in the experiments. On the other hand, short 
range interactions severely restrict the extent of the correlations, and render the dynamics 
essentially diffusive on all length scales, preventing the system from showing the collective be- 
havior reported in the experiments. We have checked this last point by performing simulations 
of a system with real walls at distances comparable to the average interparticle separation. In 
particular, with our initial conditions one obtains a sum of modified Bessel functions which 
decay exponentially fast for length scales greater than L y . As expected, after a short ballistic 
transient, we observe an essentially diffusive behavior, quite different indeed from the data 
reported in Refs.(7[ [| O. We thus conclude that long range interactions must be preserved 
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Fig. 1 - a) Snapshot of the velocity fluctuations, showing both swirls and channels, b) Trajectories 
of a fast (left) and a slow (right) particles (see text). Units given in particle radii. 



in order to account for the scale competition observed in the system. Our model is based in 
this simple consideration. 

To compute XJ n , we resort to the Ewald summation method |14[ ] , which yields the following 
expression: 



M(r) = ^H (r + d) 
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Here the function erfc(a;) is the complementary error function, and (3 is a parameter which 
controls the convergence of both the d and G sums. The reciprocal space vectors G are such 
that G • d — 2nk, where k is an integer. The terms proportional to I are the same as for the 
Coulomb potential |Q; the other terms are intrinsic to the Oseen tensor. As in the Coulomb 
case, the term G = in Eq. (|l|) yields a divergent contribution. In the electrostatic case, this 
infinite contribution cancels out after imposing an overall charge neutrality condition. In the 
sedimentation problem, the G = term cancels out after subtracting the mean sedimentation 
velocity U = (N/V) Jd 3 rHo(r) • g, and thus working with velocity fluctuations u„ = 
U„ — U M . By doing so, the particle positions are on average fixed, as in the sedimentation 
experiments in a fluidized bed. 

To follow the evolution of the trajectories and velocities of N particles, we integrate nu- 
merically the 2N coupled equations dr n /dt = J2 m M(r n — r m ) ■ g — U M , where M is given by 
Eq. ([j]), using an adaptive step-size fifth-order Runge-Kutta algorithm ]l5| . We have chosen 
a convergence parameter (3 = L/12; other values of (3 were also tested, yielding equivalent 
results. Simulations start from a configuration of N particles randomly placed on a square cell 
of size L. Since the Oseen approximation is not valid at short distances, to avoid singularities 
in the velocity field we have introduced an ad hoc very short range repulsive hard-core term 
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Fig. 2 - Integrated distributions P+ and P- for the (a) u x and (b) u z velocity fluctuations in linear-log 
scale. The solid line corresponds to an integrated Gaussian distribution. 



of the form exp[— (r — 2a) / p] , where a is the radius of the particles and p is a small parameter 
that we select equal to 0.1. 

The concentrations described by our model are severely limited by both the range of 
validity of the Oseen approximation and the available CPU time. In our simulations, therefore, 
we have considered concentrations c < 1%, and cell sizes ranging from L = 100a — 200a. We 
shall see, however, that our results for dilute concentrations already exhibit most of the salient 
features reported in the literature. Averages were made over at least 100 realizations starting 
with different random initial conditions. 

In the evolution of our model, particles build up a complex and highly fluctuating pattern of 
velocity swirls and channels, very similar to those experimentally observed j?], |). In Fig- 0a), 
we show a snapshot of a system with concentration c = 1% and cell size L — 200a. The 
number of swirls and their sense of flow (clock- or anti-clock-wise) result from the collective 
interactions, and have the constraint of zero global vorticity Y^ n V x u„ = 0, as follows from 
the symmetries of the Oseen tensor. 

At large times, the average root-mean-square velocity fluctuations (RMSVF) along the 
horizontal, u x , and vertical directions, u z , grow with the concentration c. As naively expected 
from the symmetry breaking induced by gravity, fluctuations are anisotropic. We measure 
a ratio u z /u x ~ 2.5, which seems independent of c or L. This observation agrees with the 
results reported in Ref. @j. 

Next, we have measured the probability density function (PDF) of the velocity fluctuations, 
p(u z ) and p(u x ), normalized as to have zero mean and unity standard deviation. In Fig. |^ we 
plot the integrated distribution functions P+(u) — J p(u')du' for the downward (rightward), 
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Fig. 3 - Velocity autocorrelations as a function of time. The curves shown in the inset correspond to 
the raw data, whereas in the main plot, time has been rescaled by the characteristic time r ~ c . 
(o) system I, (x) system II, (A) system III (see text). 



u > 0, velocity, and P-{— u) = f_ p(u')du' for the upward (leftward), u < 0, velocity, for 
both vertical and horizontal fluctuations. In particular, the plots correspond to values of 
c = 1% and L = 200a. We observe that the horizontal fluctuations are left-right symmetric 
and very well approximated by a Gaussian distribution (solid line in Fig. ||a)). On the other 
hand, vertical fluctuations are fairly asymmetric and apparently non-Gaussian. 

We now turn our attention to the two-times statistical properties of the velocity fluctua- 
tions. First, we consider the velocity autocorrelation function g a (t) = (u a (0)u a (t)) / ^u a (0) 2 ), 
for a = x, z, where the brackets denote an average over particles and realizations, at a fixed 
time t. In Fig. || we depict g a (t) for two different concentrations, c = 1% (system I) represented 
with (o), and c = 0.25% (system II) represented with (x), in a box of size L = 200, as well as 
c = 1% in a smaller box of size L — 100 (system III) which we plot with (A). The main plot 
represents our data as a function of the rescaled time ct; raw data are shown in the inset. For 
both g x and g Zl we observe an initial exponential decay of the correlations with a character- 
istic time proportional to c _1 . This scaling of g a at short times can be understood by means 
of a simple mean-field-like argument: Given the expression of the Oseen tensor, the velocity 
correlations can be written as (u(i)u(O)) ~ (u(0)/r(t)}, where r is the separation between any 
pair of particles. Taking a time derivative, d t (u(t)u(0)} ~ d t (u(0)/r(t)) ~ — (u(t)u(0)/r 2 ), 
where in the second step we have commuted derivative and average. A further simplifi- 
cation considers r ~ c -1 / 2 , i.e., the average separation between particles. Then, we have 
d t (u(t)u(0)} ~ — (u(t)u(0)) /c _1 , yielding an exponential relaxation with characteristic time 

T~C-\ 

After this initial decay, the x correlations of the more concentrated system I show a 
clear negative region. Curiously, this behavior resembles that of a dense liquid. Negative 
autocorrelations in a dense liquid are due to backscattering effects after collisions among 
molecules. In our system, however, negative correlations are due to the permanence of the 
particles in a velocity swirl. As argued in Q, during the course of a simulation some of the 
particles become part of velocity swirls and spend in them a considerable amount of time. 
They can be called slow particles and describe coil-like trajectories. Others (fast particles) 
spend more time inside the channels separating swirls, and their trajectories are much more 
elongated. Both channels and swirls can be observed in Fig. |l]a). In Fig. If}), we plot typical 
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trajectories of a fast and a slow particles. 

At later times the correlations of the x components oscillate around zero, whereas the z 
autocorrelations go through a second regime of much slower relaxation, and eventually become 
zero towards the end of the simulation time. This enhancement of the z autocorrelations is 
due to the very existence of channels between swirls, inside of which particles follow ballistic 
trajectories with small fluctuations. Channels are interrupted by swirls, but since these must 
be created in pairs of opposite vorticity (due to vorticity conservation) , their creation is costly 
and only a few are present in a box of small size. A long time is thus required for the particles 
initially in a channel to become part of a few velocity swirls and uncorrelate from their initial 
conditions. 

To further explore the behavior of the system at long times, we have also studied the 
mean-square displacement of the particles, R a {t) — ([r a {t) ~ r a(0)] 2 ), which is an efficient 
indicator of a possible effective diffusive behavior of the system (hydrodynamic diffusion) H . 
For the latter, we expect R a (i) — 2D a t, i.e., dR a (t)/dt = 2D a = const, where D a is an 
effective diffusion coefficient. In Fig. f| we represent the time derivative of R a (t) for the 
displacements along the x and z directions. The plateau at long times clearly indicates that 
the displacement along the x direction becomes purely diffusive right after an initial ballistic 
regime (R x (t) ~ t 2 ). The z displacement also becomes eventually diffusive, but at longer times 
scales. This final diffusive behavior is compatible with the fast decay of the tails of the PDF's 
shown in Fig. |[ We observe that D z ^> D x , hence the diffusive regime is highly anisotropic. 
At intermediate times, we observe that R z {t) can be fitted to a power-law R z (t) ~ t a with an 
exponent within the range 1 — 2. Such behavior was reported in [^), where experiments could 
not be run for long enough times as in our simulations. We expect that experiments carried 
out over longer time scales would also show the eventual diffusive behavior along the vertical 
direction. 

To sum up, we present a model for the sedimentation of non-Brownian particles in an 
unbounded fluid that incorporates long range hydrodynamic interactions and PBC's in the 
simplest Oseen approximation. This model exhibits most of the salient features of the ex- 
periments reported in Refs. jj], 0, [fl). Our findings can be understood within the picture of 
slow and fast particles: Slow particles spend most of the time within velocity swirls and con- 
tribute to the fast relaxation of the velocity correlations. Fast particles moving along velocity 
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channels have strongly correlated (quasi-ballistic) trajectories and are responsible for the slow 
relaxation component. For sufficiently long times, all particles become part of enough velocity 
swirls, and our model predicts that the system eventually relaxes to a hydrodynamic diffusive 
regime, that could be confirmed by experiments performed over longer time spans. 
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